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Abstract. Critical, or scale independent, systems are so ubiquitous, 
that gaining theoretical insights on their nature and properties has many 
direct repercussions in social and natural sciences. In this report, we start 
from the simplest possible growth model for critical systems and deduce 
constraints in their growth : the well-known preferential attachment prin- 
ciple, and, mainly, a new law of temporal scaling. We then support our 
scaling law with a number of calculations and simulations of more com- 
plex theoretical models : critical percolation, self-organized criticality 
and fractal growth. Perhaps more importantly, the scaling law is also 
observed in a number of empirical systems of quite different nature : 
prose samples, artistic and scientific productivity, citation networks, and 
the topology of the Internet. We believe that these observations pave 
the way towards a general and analytical framework for predicting the 
growth of complex systems. 
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1 Introduction 

An interesting duality of nature is that the most complex systems usually obey 
the most simple rules. For instance, fractals exhibit delicate entanglement of ge- 
ometric features which emerge from specific local rules. Similarly, the ubiquity of 
scale independent (or scale-free) organization, found throughout natural, tech- 
nological and human systems [T] , though somewhat puzzling, can be reproduced 
by simple stochastic models (e.g. critical percolation [2], self-organized criticality 
[3 , diffusion-limited aggregation [4] and preferential attachment [5]). Our project 
aims to uncover the secret ingredients shared by these stochastic models and the 
empirical data they reproduce, in the hope of unifying our understanding of the 
growth of complex critical systems. 

More precisely, the present report deals with the temporal evolution of one 
particular feature of these critical systems: the scale independent distribution of 
a given quantity K among the N elements of the system. 

Mathematically, two very general hypotheses will be used to describe the 
systems under study: 
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1. the asymptotic proportion nk{t — > oo) of elements of size k follows a power 
law: — > oo) oc k 1 (with 7 > 1 for normalization); 

2. relative to other elements, an element will grow at a rate G(k) dictated 
mostly by its size k and not by some hidden variable. 

As a consequence of these two hypotheses, we can state that G(k) oc A; at least for 
k 3> 1 (see Appendix A). The second hypothesis also implies that elements of a 
given size k are indiscernible from each other. As a general notation, we will use 
Nk(t) for the average number of elements of size k at time t, such that N(t) = 
J2 k Nk(t) is the average total population at time t, and n k {t) — Nk(t)/N(t), the 
studied distribution. We can now state that the growth of the system follows a 
simple set of equations: 



The first p(t) is the probability of a birth event (which affect size k — 1 only), 
and the last term is the probability of a growth event affecting an element of size 
k (either by bringing an element of size k — 1 to k or by moving an element of 
size k to k + 1). All variables are discrete as size is numbered in units of growth 
events, just as time t corresponds to the total number of events (birth or growth) 
since the birth of the system. 

Our hypothesis # 2 provides a general form for the growth junction, G{k) oc k, 
which is typically referred to as preferential attachment. To completely define 
Eq. 0, we must find a similarly general form for the birth function pit). This 
function yields the probability that the t-th event in the system's history corre- 
sponds to the birth of a new element. The search for a functional form, analytic 
or heuristic, of pit) will occupy a major part of our contribution. Note that Eq. 

does not consider the potential death of elements, these events are considered 
indiscernible from the state of "not growing ever again" and are thus averaged 
within the growth function for simplicity. 

2 An approximate law 

We study the probability pit) that the t-th. event in the system marks the birth 
of a new element, as opposed to a growth event. Let us consider a system which 
follows the simplest growth rules, i.e., apart from our two hypothesis, let us 
assume that the birth rate is given by AiN{t) while the growth rate of an 
element of size k is BiGik). The system is conservative (i.e. no death events) 
such that we are interested in how pit) falls to zero as the population goes to 
infinity. We can directly write p[t) as the ratio of birth rate to the sum of birth 
and growth rates: 



N h (t + 1) = N k (t)+p(t)5 kl + (l- p(t)) 



[G jk - 1) Nk_.it) - G jk) N k jt)} 
Ek>G(k')N k >(t) 



(1) 



p(t) 



AiNit) 



(2) 



AiN(t) + E 




SiAT(i)n fe (t)G(fc)dA; 
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As rik(t)G(k) is a monotonous and slowly decreasing function of k, we can ap- 
proximate the finite sum by the equivalent integral and relegate the error made 
for very small k as part of the integration constamQ Equation ^ becomes 

Pit) ~ Mm 



A 1 N{t) + ^2T {t) B 1 N{t)n k {t)G{k)dk 
* 1 = 1 ( 3 ) 

i + B 2 J^ (t) k-^kdk i^k max (t) 2 -y + C 

where all capital case letters are irrelevant constants. To further simplify this last 
expression, one must determinate the behaviour of the largest element k max (t) 
at time t. In Appendix B, we show that fcmax(i) — Pot 1 ~ p ^ t ' + pi, which yields 

P ^ ~ C + D^-pW + £) 2 -t ' (4) 

As pit) falls towards zero, we see that it reaches an asymptotic temporal 
scaling with (t + E) with exponent cx (2 — 7) . Note that assuming normalization 
of the target distribution in k^ 1 implies 7 > 1, thereby allowing the temporal 
dependence of p(t) to scale with either positive or negative exponents. However, 
because of the challenge inherent to accurately evaluating the slope of a power- 
law [B], especially if it changes as the system evolves, we will leave this exponent 
adjustable. 

Furthermore, as the only term function of t in Eq. Q is related to the total 
number of growth events, it is expected to be much larger than C. Hence, it is 
useful to approximate the global behaviour of p(t) with the functional form: 

p(t) = a{t + r) a + b (5) 

where r delays the temporal scaling in a towards the limit b. We have here 
assumed that dissipative systems (with death events) will follow a qualitatively 
similar dynamics as p(t) goes to its asymptotic value b. This hypothesis will be 
tested on various systems. 



3 Critical percolation 

Let us now consider a two-dimensional square lattic^] where at each step (i —> 
i + 1) one of L sites is randomly chosen for a percolation trial. With probability 
u, the site is occupied and the system clock is increased by one (t — > t+ 1). Sites 
can only be chosen once for this trial. 

1 The error due to this continuous approximation goes as k~~ 2 and is therefore signif- 
icant for small k solely. This implies that the error due to the approximation of the 
sum by the integral is independent of the growth of fc ma3E (t). 

2 Each site of a square lattice has exactly four neighbours. We consider the infinite 
size limit to neglect finite size effects. 
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We study the behaviour of connected clusters. What is the probability p(t) 
that the t-th occupied site creates a new cluster of size one (none of his four neigh- 
bours is occupied yet)? This probability defines the ratio between the growth of 
the cluster population to the number of occupied sites. We expect p(t) to follow 
Eq. ([5]) , since it is known that percolation leads to a critical state where cluster 
sizes follow a power-law distribution (hypothesis 1) and it is simple to imagine 
how a larger cluster has more growth potential than a smaller one because it has 
more neighbouring sites (hypothesis 2^ 

The t-th occupied site will mark the birth of a new cluster if none of his 
four neighbours were among the (t — 1) first occupied sites. We can then directly 
write: 

- ri f 1 - L=j) • ( g ) 

Rewriting p(t) as 



one can see that 

4 

p{t) = n L ~ t L + _ 1 ~ i f ° r t > 4 ( § ) 

i=i 3 

For L » t » 1, 

Equation ^ agrees with Eq. ^ using a = L~*, r = —L, a = 4 and 6 = 0, see 
Fig. [T] It is not surprising to see a > as critical percolation in two dimensions 
leads to 7 = 187/91 > 2 such that 2 - 7 < in Eq. Q. 

It is important to note that Eq. ^ docs not depend on the percolation 
probability u. Under this form, the ability of this system to converge towards 
its critical state depends on the number of sites occupied, i.e. time t. Noting 
that t = uL, the critical time t c = u c L, corresponding to the critical point in it, 
could perhaps be calculated through a self-consistent argument on the number of 
occupied sites required by the scale-free distribution of cluster size in the critical 
state. This is however, not our current concern. 



4 More complex processes 

This section presents numerical studies of growth processes for which analytical 
solutions of the birth function pit) are unavailable. 

3 In fact, it can be shown that in the limit of size k ^> 1, the perimeter of a cluster of 
size k scales as k [7]. 



Critical Growth 5 




200000 400000 10° 10 1 10 2 10 3 10 4 10 5 

t t 



(a) (b) 

Fig. 1. Percolation on a lOOOxlOOO square lattice at phase transition 
(u c = 0.5927...). (a) The evolution of the probability p(t) that the t-ih occupied 
site results in the creation of a new cluster (in semi-log plot), (b) Log-log plot of 
the complementary probability 1 — p(t) to highlight the initial temporal scaling. The 
solution corresponds to Eq. Q9p . 



4.1 Self-organized criticality: sandpile model 

In constrast to systems, like percolation, which are critical at a phase transition, 
other systems are known to always converge to their critical state, without the 
need to adjust any parameters. This asymptotic criticality is often referred to as 
self- organized criticality (SOC). The first and simplest example of SOC is the 
Bak-Tang-Wiesenfeld model [8], a sandpile system where grains are dropped on 
a surface and grains Zi of a given site i topple on the four nearest neighbours 
z nn if Zi reaches a threshold z* = 4. The algorithm is 

1. Initialisation. Prepare the system in a stable configuration: we choose Zi — V i. 

2. Drive. Add a grain at random site i. 

Zi — S> Zi + 1 . 

3. Relaxation. If z% > z* , relax site i and increment its 4 nearest-neighbours (nn). 

Zi — > Zi — 4 , 

Znn ' Z n n 

Continue relaxing sites until z% < z* for all i. 

4. Iteration. Return to 2. 

In addition to the self-organized property of this model (hypothesis 1), larger 
potential avalanches have more neighbouring sites resulting in a bigger growth 
potential (hypothesis 2). Hence, the probability that the i-th site marks the 
creation of a new potential avalanche is expected to follow Eq. Moreover, we 
can also expect the constant b to be non-zero because of the dissipative nature 
of the model: avalanches occur and thus stop growing, new potential avalanches 
will consequently continuously appear. Results of simulations of the BTW model 
are presented on Fig. [2j 
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Fig. 2. BTW sandpiles on a 64x64 square lattice with z* = 4. (a) The evolution 
of the probability p(t) that the t-th site to reach a height of z* — 1 results in the creation 
of a new potential avalanche (in log-log plot), (b) Plot of the complementary probability 
1 - p(t) for growth events. The fit uses Eq. (5} with b = 0.3090, a = -3.5, r = 3000 
and a fixed by p(l) = 1 (i.e. a = (1 - 6)(1 + i") - "). 



Fig. 3. Fractal structure created through diffusion-limited aggregation. Re- 
sult of one realization of diffusion-limited aggregation on a square lattice (512x512) 
with an initial circular structure (radius r = 82) at its center. 



4.2 Fractal growth: diffusion-limited aggregation 

Diffusion-limited aggregation (DLA) is a stochastic process where particles, un- 
dergoing Brownian motion, cluster together to form aggregates [J| (see Fig. [3]). 
The process was developed as a model of aggregation in systems where diffusion 
is the primary means of transport and where forces are negligible. The resulting 
aggregates have been shown to follow a fractal structure where density corre- 
lations fall off with distance as a fractional power law (dimension ~ 1.7). Real 
examples of DLA-like aggregation can be observed in systems such as electrode- 
position, mineral deposits and dielectric breakdown. 

Using results of Hinrichsen et al. [S], it is easy to show that the length of 
branches in the Brownian trees produced by DLA follows a power-law distri- 
bution. Here, branches are differentiated according to their order: the trunk is 
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Fig. 4. Growth of a DLA fractal aggregate on a 750x750 square lattice 
with an initial circular structure (r = 82) at its center. The evolution of the 
probabilities p(t) (red) and 1 — p(t) (blue) that the t-th grain creates a new branch 
or that it continues an existing branch, respectively, in numerical simulations (log-log 
plot). The fit uses Eq. (5) with b = 0.23, a = -7.5, r = 720 and a fixed by p(l) = 1. 



the zeroth-order branch and branches stemming from a branch of order n are of 
order n + 1. Results demonstrate that the length L(n) and the number N(n) of 
branches of order n are exponential functions of n: 

L(n) w d? and N(n) w d% with d x ~ 0.36 and d 2 ~ 5.2. (10) 

The number AT(L) of branches of length i is found using n ~ logL/logcii: 

N(L)^L~~ ( with 7 = -logd 2 /logdi. (11) 

The distribution of the lengths of the branches is scale-free most likely because 
longer branches gain an advantage by overshadowing smaller branches, i.e. the 
higher a branch reaches, the more likely it is that a random trajectory crosses 
its path (hypothesis 2). We will follow the probability p(t) that a grain which 
reaches the aggregate marks the beginning of a new branch by attaching to the 
body of a branch and not to its tip. Figure [4] illustrates how this property evolves 
according to Eq. ^ The system is not, strictly speaking, conservative, because 
branches can lose the ability to grow once it becomes impossible for random 
trajectories to reach their tip. 



5 Empirical systems 

5.1 Word occurrences in prose samples 

This subsection is concerned with the reading of prose samples. It is well-known 
that the word frequency distributions of written texts are approximately scale- 
free, i.e. the number of words which appear k times in a text falls roughly as 
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Fig. 5. Prose samples from William Shakespeare's work. The evolution of the 
probability p(t) that the t-th written word is a new word in the vocabulary used within 
this sample (log-log plot). The fit uses Eq. |5| with b — 0, a — —0.21, r = 22 and a 
fixed by p(l) = 1. 

fc~ 7 [TO]. We shall study here the probability p(t) that the t-th word is a word 
that has yet to appear in the text. By its empirical nature, this system is a 
bit more complicated than the previously considered theoretical systems. Some 
words can be expected to follow hypothesis 2 (e.g. adjectives and nouns which 
refer to the main subject or setting), while others (e.g. determinants) will appear 
more frequently as syntaxic constraints. Hence, the system can be expected to 
behave as a hybrid between the critical systems discussed above and the mere 
random sampling of a scale-free system. 

Using all samples of 1000 consecutive words in the writings of William Shake- 
speare (around 9 x 10 5 words) [IT] , we can take the mean value of p(t) across 
samples for each t £ [1,1000]. Results of this simple experiment are presented 
on Fig. [5] 

5.2 Scientific and artistic productivity 

Most scale-free empirical data pertain to unique systems, for which it is impos- 
sible to obtain the mean function p(t) as we did for prose samples. Instead, the 
history of these systems is a binary sequence, i.e. p{t) is 1 if the t-th event is a 
birth event and if it is a growth event. To obtain a continuous p(t), we simply 
use a running average procedure with windows of At. See Fig. [6] for the result 
of this method on the datasets of authorship of scientific articles in the arXiv 
database [12] and of castings in the Internet Movie Database [13j . 

5.3 Complex networks 

This last section highlights the relation between the temporal scaling presented in 
this report and the so-called densification of complex networks. This last concept 
refers to the evolution of the ratio of links to nodes in connected systems. In the 
case of scale-free (or scale independent) networks, this densification was observed 
to behave as a power-law relation between the number of nodes and the number 
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Fig. 6. Scientific and artistic productivity, (a) The evolution of the probability 
p(t) that the t-th name to appear as an author of an article on the arXiv database 
(obtained through running average with At — 10000) is a name which appears for the 
first time. The fit uses Eq. (5} with a = 97, b = 0.235, r = 7000 and a = -0.55. (b) 
The evolution of the probability 1 — p(t) that the t-th actor cast for a role will not be a 
new actor (running average with At — 25000). The fit of our law (black) uses a = —11, 
r = 10 6 , a = —0.18 and 6 = 1. The dotted line is the asymptotic power-law behaviour 
of both the data and our fit. 



of links [14] . Based on our previous results, we can conjecture a more precise 
relation. 

In analogy to our law, the number M of links would be directly proportional 
to the total number of events (or more precisely time M — t/2 as one link involves 
two nodes) while the number of nodes is directly related to the total population 
N(t). As these databases do not consider node removal (death events), we have 
N(t) = J Q p(t')dt' and b — 0. Hence we expect the numbers of nodes and links 
to be related through the following expression 



N(M) 



{2M + t) l 



r a+l 



(12) 



Eq. (12 1 can be rewritten as 
N(M) ~ — 



(1 + 2M/t) c 



> + 1 v ' ' a + 1 

to show that the relation is initially linear, i.e. when t and tf «t, 



(13) 



iV(M<r) 



(IT 
O 



a+1 
1 



x 2M ^ 
l+(a+l) +0 



M' 



2aT a M. (14) 



Equation ( 12 ) thus predict an initially linear densification leading into a power- 



law relation. This behaviour is tested on two network databases: the topology of 
Autonomous Systems (AS) of the Internet in interval of 785 days from November 
8 1997 to January 2 2000 and the network of citations between U.S. patents as 
tallied by the National Bureau of Economic Research from 1963 to 1999 pL4j. 
The results are presented on Fig. [7] Of the four systems considered in [14] , these 
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two were chosen to highlight two very different scenarios. On the one hand, the 
Internet can be reproduced by multiple pairs of a and r parameters as long as 
r«J, since the system appears to have reached steady power-law behaviour. 
On the other hand, the patent citation networks do not fit with the power-law 
hypothesis as the system is transiting from a linear to a sub-linear power-law 
regime as t ~ r. This last scenario, while very different from a simple power-law 
growth as previously proposed, corresponds perfectly to the predictions of our 
theory. 




Fig. 7. Densiflcation of complex networks in log-log scale. Densification, i.e. 
the relation between the number of nodes and the number of links, in two connected 
systems, (a) The Internet at the level of autonomous systems, reproduced by Eq. ( 12 1 
with a ~ 1, r — 0, a = —0.16 and b = 0. (b) The citation network of U.S. patents 
between 1963 to 1999, reproduced by Eq. (12} with a = 5793, r = 875000, a = -0.67 
and 6 = 0. The dotted line is the power-law relation claimed in [14j . 



6 Conclusion and perspectives 

In conclusion, we have provided theoretical arguments and empirical justifica- 
tions that a system which converges toward a scale-free (critical) organization 
through constant growth rules will feature delayed temporal scaling in the fre- 
quency of birth events. We have introduced an approximation for the proba- 
bility p(t) that the i-th event is a birth event of the generic form of Eq. 
p(t) = a{t + T) a + b. 

This approximation was shown to be correct in various critical models: per- 
colation, self-organized criticality and fractal growth; and in empirical systems: 
prose samples, the movie industry, the Internet and the scientific literature. The 
next step to build on this theory is to add stronger analytical foundations and 
to use it as a predictive tool. Now that Eq. ([I]) is complete with functional forms 
for both the growth and birth functions, it can be used to replicate systems 
for which temporal data is unavailable. From a single snapshot of a system's 
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present state, we could hope to provide an educated guess for a system's past 
and perhaps even predict its future evolution. 



Appendix A: Growth function 

This first Appendix concerns the proof that a scale-free distribution for Nk(t) 
implies an asymptotically linear preferential attachment ; in other words that 
G(k —> oo) oc k. Let us redefine the growth rate G(k) as the average fraction of 
elements of size k which will grow during the upcoming time step. We write: 

G(k)=j— ) [Niit + V-Niit)} (15) 

where the sum yields the number of elements that left the compartment dur- 
ing this time step and the normalization gives the desired fraction. As mentioned 
in the main text, the sum can be replaced by an integral when dealing with large 
values of k. Hence, it is easily obtained that for Nk(t) = A(t)k~' y 

G(k)~ A{t + ]}~ A(t) k ocfe forfc^oo. (16) 
v ; A(t) 7-1 v ' 



Appendix B: Growth of maximal size 



This second Appendix proves that, assuming a linear growth function G(k) ~ k 
and a slowly varying birth function p(t) , the maximal element size /c max (t) present 
in the system at time t scales as t~@ [15] . In fact, because Eq. jl]) is deterministic, 
the first element is certain to follow fc max (i). Moreover, the chosen G(k) implies 
that the normalization of growth rates at time t will always be approximately 
given by t. Denoting ki(t) the size of an element i at time £, we can write the 
result of a single step of Eq. as: 



h{t+i) 



p (t) + (i- P (t)) 



t-ki(t) 



t 



ki(t)- 



(i-K*)) 



Simplifying the equation yields 

ki(t+l) = 



l-p(t) 



t 



ki(t) 



(hit) + 1) . (17) 



(18) 



which fixes the derivative in the limit of large t: 



(19) 



Assuming that lim^oo N(t) — > oo, requires p(t) to either converge toward a 
constant, or fall slower than t . The two options imply that the solution to Eq. 



(19) can be approximated by the following Ansatz: 
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ki(t) a at 1 -"® + ca (20) 

which is a general solution for any k(t), where C\ and C2 depend on p(t) and 
initial conditions. 



Appendix C: Summary of produced fits 



process 




a 




r 


a 


b 


percolation 




L" a 




— L 


4 





SOC (BAK sandpiles) 


(1 


-6)/(r- 




~ L 


-3.5 


0.31 


Fractal growth (DLA) 


(1 


-6)/(r- 




~ L 


-7.5 


0.23 


Word occurrences 




(1+r)- 


a 


22 


-0.21 





arXiv authorship 




97 




7000 


-0.55 


0.235 


IMDb castings 




-11 




10 6 


-0.18 


1 


Internet structure 




1 







-0.16 





patents citations 




5793 




875000 


-0.67 






Table 1. Summary of parameters used in Eq. ([5J) to fit diverse processes. 

Note: L is the lattice size used in corresponding simulations. 
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